//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// Compare interpolation to direct calculation

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"

#define RADIUS 0.15
#define WR 0.035
#define SPACING 0.01
#define I 50000
#define Q 3e-7

int main (int argc, char *argv[])
{
  FILE *fp = fopen ("ipoly.csv", "w");
  if ( !fp )
    {
      fprintf (stderr, "Unabel to open ipoly.csv for writing!\n");
      exit(-1);
    }

  fprintf (fp, "MEMORY[MB],TICKS\n");

  for ( size_t i = 1; i <= 10; i++ )
    {
      {
      CoilGen cgen;
      cgen.radius = RADIUS;
      cgen.wr = WR;
      cgen.set_space (SPACING);
      cgen.do_efield = true;
      cgen.shape = CoilGen::ROUND;
      cgen.geom_radial = 6 * RADIUS;
      cgen.geom_axial = 6 * RADIUS;
      cgen.memory_usage = i * 10;
      cgen.coil_matrix_n = 1;

      printf ("Running test with %lu MB\n", (unsigned long)cgen.memory_usage);
      cgen.generate ("ipoly.cgen", false, false);
      }

      {
      CoilGen cgen;
      cgen.radius = RADIUS;
      cgen.wr = WR;
      cgen.set_space (SPACING);
      cgen.do_cubic = true;
      cgen.do_efield = true;
      cgen.shape = CoilGen::ROUND;
      cgen.geom_radial = 6 * RADIUS;
      cgen.geom_axial = 6 * RADIUS;
      cgen.memory_usage = i * 10;
      cgen.coil_matrix_n = 1;

      printf ("Running test with %lu MB\n", (unsigned long)cgen.memory_usage);
      cgen.generate ("ipoly2.cgen", false, false);
      }

      CoilData cd;
      if ( !cd.load ("ipoly.cgen") )
        {
          fprintf (stderr, "Unable to load generated ipoly.cgen\n");
          exit(-1);
        }

      CoilData cd2;
      if ( !cd2.load ("ipoly2.cgen") )
        {
          fprintf (stderr, "Unable to load generated ipoly2.cgen\n");
          exit(-1);
        }

      Statics dcs, nns, lis, cus;
      make_polywell_cube (dcs, RADIUS, WR, SPACING, I, Q, ShapedCoilFactory());
      make_polywell_cube (nns, RADIUS, WR, SPACING, I, Q, CoilDataFactory<StaticCINN>(cd));
      make_polywell_cube (lis, RADIUS, WR, SPACING, I, Q, CoilDataFactory<StaticCIL>(cd));
      make_polywell_cube (cus, RADIUS, WR, SPACING, I, Q, CoilDataFactory<StaticCICU>(cd2));

      ElectroDynamics dcdyn(dcs);
      ElectroDynamics nndyn(nns);
      ElectroDynamics lidyn(lis);
      ElectroDynamics cudyn(cus);

      vect3d pos(-2*RADIUS, 0.015, 0.015);
      vect3d dir(1, 0, 0);
      prec_t speed = 0;

      vect3d bf, ef;
      dcs.getFields (pos, bf, ef);
      bf.print("dcs.bf");
      ef.print("dcs.ef");
      nns.getFields (pos, bf, ef);
      bf.print("nns.bf");
      ef.print("nns.ef");
      lis.getFields (pos, bf, ef);
      bf.print("lis.bf");
      ef.print("lis.ef");
      cus.getFields (pos, bf, ef);
      bf.print("cus.bf");
      ef.print("cus.ef");

      dcdyn.inject_electron (pos, dir, speed);
      nndyn.inject_electron (pos, dir, speed);
      lidyn.inject_electron (pos, dir, speed);
      cudyn.inject_electron (pos, dir, speed);

#define SCREEN_SIZE 380

      Screen dcscr(SCREEN_SIZE, SCREEN_SIZE, 4.5*RADIUS, 4.5*RADIUS);
      Screen nnscr(SCREEN_SIZE, SCREEN_SIZE, 4.5*RADIUS, 4.5*RADIUS);
      Screen liscr(SCREEN_SIZE, SCREEN_SIZE, 4.5*RADIUS, 4.5*RADIUS);
      Screen cuscr(SCREEN_SIZE, SCREEN_SIZE, 4.5*RADIUS, 4.5*RADIUS);
      dcscr.draw_axes ();
      nnscr.draw_axes ();
      liscr.draw_axes ();
      cuscr.draw_axes ();
      dcs.draw (dcscr);
      nns.draw (nnscr);
      lis.draw (liscr);
      cus.draw (cuscr);

      prec_t maxdev = 0.01 * 0.01;

      prec_t cu_max_erre, cu_max_errb, cu_cum_erre, cu_cum_errb;
      prec_t li_max_erre, li_max_errb, li_cum_erre, li_cum_errb;
      cu_max_erre = cu_max_errb = cu_cum_erre = cu_cum_errb = 0;
      li_max_erre = li_max_errb = li_cum_erre = li_cum_errb = 0;
      size_t n = 0;
      while (1)
        {
          dcdyn.step();
          nndyn.step();
          lidyn.step();
          cudyn.step();

          dcscr.set_pixel (dcdyn.get_pos(0), Screen::RED);
          nnscr.set_pixel (nndyn.get_pos(0), Screen::RED);
          liscr.set_pixel (lidyn.get_pos(0), Screen::RED);
          cuscr.set_pixel (cudyn.get_pos(0), Screen::RED);
          if ( (dcdyn.get_pos(0) - cudyn.get_pos(0)).magnitude() > maxdev )
              break;

          vect3d bf1, ef1, bf2, ef2;
          cus.getFields (dcdyn.get_pos(0), bf1, ef1);
          dcs.getFields (dcdyn.get_pos(0), bf2, ef2);
          prec_t cu_bfdiff = (bf1-bf2).length();
          prec_t cu_efdiff = (ef1-ef2).length();
          prec_t cu_errb = cu_bfdiff / bf2.length();
          prec_t cu_erre = cu_efdiff / ef2.length();
          cu_cum_errb += cu_errb;
          cu_cum_erre += cu_erre;
          if ( cu_max_errb < cu_errb ) cu_max_errb = cu_errb;
          if ( cu_max_erre < cu_erre ) cu_max_erre = cu_erre;

          cus.getFields (dcdyn.get_pos(0), bf1, ef1);
          lis.getFields (dcdyn.get_pos(0), bf2, ef2);
          prec_t li_bfdiff = (bf1-bf2).length();
          prec_t li_efdiff = (ef1-ef2).length();
          prec_t li_errb = li_bfdiff / bf2.length();
          prec_t li_erre = li_efdiff / ef2.length();
          li_cum_errb += li_errb;
          li_cum_erre += li_erre;
          if ( li_max_errb < li_errb ) li_max_errb = li_errb;
          if ( li_max_erre < li_erre ) li_max_erre = li_erre;

#if 0
          char buf[128];
          printf ("%s: %g (%.8f%%), %g (%.8f%%)\n",
              dcdyn.get_pos(0).sprint(buf), cu_bfdiff, 100 * cu_errb, cu_efdiff, 100 * cu_erre);
#endif

          n++;
        }
      printf (" cubic field error: max_errb=%.8f%% (avg %.8f%%), max_erre=%.8f%% (avg %.8f%%)\n",
          prec2double(100*cu_max_errb), prec2double(100*cu_cum_errb / n),
          prec2double(100*cu_max_erre), prec2double(100*cu_cum_erre / n));
      printf ("linear field error: max_errb=%.8f%% (avg %.8f%%), max_erre=%.8f%% (avg %.8f%%)\n",
          prec2double(100*li_max_errb), prec2double(100*li_cum_errb / n),
          prec2double(100*li_max_erre), prec2double(100*li_cum_erre / n));

#if 0
      vect3d p(-0.000812064, 0.00210751, 0.000679883);
      vect3d b1, b2, b3;
      vect3d e1, e2, e3;
      dcs.getFields (p, b1, e1);
      cus.getFields (p, b2, e2);
      lis.getFields (p, b3, e3);
      b1.print("direct.b");
      b2.print("cubic.b");
      b3.print("linear.b");
      e1.print("direct.e");
      e2.print("cubic.e");
      e3.print("linear.e");
      printf ("b.error = %.5f%%\n", 100.0 * (b1 - b2).length() / b1.length());
      break;
#endif

      char name[128];
      sprintf (name, "ipoly_dc_%04lu.bmp", (unsigned long)(i * 10)); dcscr.write (name);
      sprintf (name, "ipoly_nn_%04lu.bmp", (unsigned long)(i * 10)); nnscr.write (name);
      sprintf (name, "ipoly_li_%04lu.bmp", (unsigned long)(i * 10)); liscr.write (name);
      sprintf (name, "ipoly_cu_%04lu.bmp", (unsigned long)(i * 10)); cuscr.write (name);

      printf ("Simulation ended in %lu ticks\n", (unsigned long)dcdyn.get_ticks());

      fprintf (fp, "%lu,%lu\n", (unsigned long)(i * 10), (unsigned long)dcdyn.get_ticks());
      fflush (fp);
    }

  fclose (fp);

  return 0;
}

